%% Optimize JRD model
clc
clear
rng('shuffle')

global newTestProbe
global obsRespErr
global refDirectionFull

% Prepare data
% Individual fitting
Exp2matrix = readmatrix('Experiment 2 (LayM).xlsx');
Exp2objs = convertCharsToStrings(readcell('Experiment 2 (LayM) objects.xlsx'));

% Inidividual fitting
for i=1:3
    idx1 = find(Exp2objs(:,i) == "JAR");
    idx2 = find(Exp2objs(:,i) == "PHONE");
    idx3 = find(Exp2objs(:,i) == "CLOCK");
    idx4 = find(Exp2objs(:,i) == "BOOK");
    idx5 = find(Exp2objs(:,i) == "BRUSH");
    idx6 = find(Exp2objs(:,i) == "CAP");
    idx7 = find(Exp2objs(:,i) == "SHOE");
    idx8 = find(Exp2objs(:,i) == "CAN");
    idx9 = find(Exp2objs(:,i) == "BANANA");
    Exp2matrix(idx1,i+5) = 1;
    Exp2matrix(idx2,i+5) = 2;
    Exp2matrix(idx3,i+5) = 3;
    Exp2matrix(idx4,i+5) = 4;
    Exp2matrix(idx5,i+5) = 5;
    Exp2matrix(idx6,i+5) = 6;
    Exp2matrix(idx7,i+5) = 7;
    Exp2matrix(idx8,i+5) = 8;
    Exp2matrix(idx9,i+5) = 9;
end

% Initial parameter values
x0 = [.4 4];

% Use this for global search
opts = optimset('Algorithm','interior-point');
problem = createOptimProblem('fmincon','options',opts,'objective',@JRDLogGrowthModelFullRefFixL, ... 
    'x0',x0,'lb',[0,0]);
gs = GlobalSearch;

for p=1:1
    
refDirectionFull = 1;

for i=1:14
    newTestProbe = Exp2matrix(1+(128*(i-1)):(128*i),6:8);
    obsRespErr = Exp2matrix(1+(128*(i-1)):(128*i),13);
    optFull(i,:) = run(gs,problem);
end

% Median and mean parameters for full reference directions
med_optFull_params(1,1) = median(optFull(:,1));
med_optFull_params(1,2) = median(optFull(:,2));
mean_optFull_params(1,1) = mean(optFull(:,1));
mean_optFull_params(1,2) = mean(optFull(:,2));


refDirectionFull = 0;

for i=15:28
    newTestProbe = Exp2matrix(1+(128*(i-1)):(128*i),6:8);
    obsRespErr = Exp2matrix(1+(128*(i-1)):(128*i),13);
    optHalf(i-14,:) = run(gs,problem);
end

% Median and mean parameters for 0 and 180 reference direction
med_optHalf_params(1,1) = median(optHalf(:,1));
med_optHalf_params(1,2) = median(optHalf(:,2));
mean_optHalf_params(1,1) = mean(optHalf(:,1));
mean_optHalf_params(1,2) = mean(optHalf(:,2));

for i=1:10000
    indx = randi([1 14]);
    bs_optFull(i,:) = optFull(indx,:);
    bs_optHalf(i,:) = optHalf(indx,:);
end

% Median and mean parameters for full reference directions (bootstrap)
bsmed_optFull_params{p}(1,1) = median(bs_optFull(:,1));
bsmed_optFull_params{p}(1,2) = median(bs_optFull(:,2));
bsmean_optFull_params{p}(1,1) = mean(bs_optFull(:,1));
bsmean_optFull_params{p}(1,2) = mean(bs_optFull(:,2));

% Median and mean parameters for 0 and 180 reference direction (bootstrap)
bsmed_optHalf_params{p}(1,1) = median(bs_optHalf(:,1));
bsmed_optHalf_params{p}(1,2) = median(bs_optHalf(:,2));
bsmean_optHalf_params{p}(1,1) = mean(bs_optHalf(:,1));
bsmean_optHalf_params{p}(1,2) = mean(bs_optHalf(:,2));

end
